Clustering Notes

Clustering Notes

The following notes follow along this blog post. Clustering methods all use a central idea of a distance or metric defined by a function \(d(i,j)\). It is often the Euclidean distance.

DBSCAN

Concepts

Consider a set of points in some space to be clustered. Let \(X = \{x_1, ..., x_n\}\), a set of \(n\) objects, \(\epsilon\) be a parameter specifying a radius of a neighborhood with respect to some point. Points in DBSCAN fall into 3 categories: - a point \(p\) is a core point if there are at least a minimum number of points within \(\epsilon\) of it, including itself - a point \(q\) is directly reachable from point \(p\) if \(q\) is within distance \(\epsilon\), but only if \(p\) is a core point - a point \(q\) is reachable form \(p\) if there is a path \(p_1, ..., p_n\) with \(p_1 = p\) and \(p_n = q\) where each \(p_{i+1}\) is directly reachable from \(p_i\). This implies that the initial point and all points on the path must be core points, with the potential exception of \(q\) - all points not reachable from any other point are outliers or noise points

Reachability is not a symmetric relation, so another concept of connection is required for the extent of clusters: - two points are density-connected if there is a point \(o\) such that both \(p\) and \(q\) are reachable from \(o\)

Clusters satisfy two properties: - all points within a cluster are mutually density-connected - if a point is density-reachable from some point of the cluster, it is part of the cluster as well

Steps

0. Choose minimum points \(m_{pts}\) threshold and \(\epsilon\)

1. Find the points in the \(\epsilon\) neighborhood of every point, and identify the core points with more than \(m_{pts}\) neighbors

2. Find the connected components of core points on the neighbor graph, ignoring all non-core points

3. Assign each non-core point to a nearby cluster if the cluster is an \(\epsilon\) neighbor, otherwise assign it to noise

Optimization Criterion

For any possible clustering \(C = \{C_1, ..., C_t\}\) out of the set of all clusterings \(\mathcal{C}\), DBSCAN minimizes the number of clusters such that every pair of points, \(p, q\) in a cluster is density-reachable.

\[\underset{C \subset\mathcal{C}, d_{db}(p, q) \le \epsilon \forall p, q \in C_i \forall C_i \in C}\min |C|\]

HDBSCAN

This is a hierarchical extension of DBSCAN. For a chosen \(m_{pts}\)

Concepts

Core Distance \(d_{core}(x_p)\)

applies to object \(x_p \in X\), with respect to \(m_{pts}\), is defined as the distance from \(x_p\) to its \(m_{pts}\)-nearest neighbor, including itself

\(\epsilon\)-Core Object

Any object \(x_p \in X\) is where the following holds: \[d_{core}(x_p) \le \epsilon\]

Mutual Reachability Distance \(d_{mreach}(x_p, x_q)\)

is defined as \[d_{mreach}(x_p, x_q) \equiv\max\{d_{core}(x_p), d_{core}(x_q), d(x_p, x_q) \}\]

Mutual Reachability Graph \(G_{m_{pts}}\)

The complete graph in which the objects of \(X\) are vertices and the weight of each edge is the mutual reachability distance (wrt \(m_{pts}\)) between the respective pair of objects.

Subgraphs of \(G_{m_{pts}}\) and Hierarchies

If we take a subgraph \(G_{m_{pts}, \epsilon}\) obtained by removing all edges from \(G_{m_{pts}}\) with weights greater than \(\epsilon\), then the connected components are the \(\epsilon\)-core objects in the DBSCAN clusters, and remaining objects noise. Thus, DBSCAN partitions for \(\epsilon \in [0, \infty)\) can be produced in a hierarchical way: at 0 we have the finest clusters, or leaves of a dendrogram, while at higher values, we approach the root of the dendrogram.

Steps

Core Algorithm

1. Compute core distance wrt \(m_{pts}\) for all data objects in \(X\)
2. Compute minimum spanning tree (MST) of \(G_{m_{pts}}\), the mutual reachability graph
3. Extend the MST to obtain \(MST_{ext}\) by adding for each vertex, a “self edge” with the core distance of the corresponding object as weight
4. Extract the HDBSCAN hierarchy as a dendrogram from \(MST_{ext}\)
   • 4.1 For the root of the tree, assign all objects the same label
   • 4.2 Iteratively remove all edges from the \(MST_{ext}\) in decreasing order of weight, with ties being removed simultaenously
     • 4.2.1 Before each removal, set the dendrogram scale value of the current hierarchical level as the weight of the edge(s) to be removed
     • 4.2.2 After each removal, assign labels to connected component(s) that contain the end vertices of removed edges, a new cluster label to a component if it still has at least one edge, else noise

Customizable Cluster Processing Final Step

When reaching a lower \(\epsilon\) level in the tree, we must consider what happens to the clusters. We can define a separate minimum cluster size which can be applied to the following: - if all cluster’s subcomponents are spurious, it disappears - if only one of cluster’s subcomponents is not spurious, keep original label, ie cluster shrinkage - if two or more of cluster’s subcomponents are not spurious, “true” cluster split

Optimal Clustering

Cluster Stability

When increasing the density threshold, we can see prominent clusters remain and shrink or split while others disappear. A cluster can be seen as a set of points that whose density \(f(x)\), exceed a threshold, ie. \(\lambda = 1/{\epsilon}\). For a density contour cluster \(C_i\) that appears at density level \(\lambda_{min}(C_i)\), we can formalize its stability by defining its excess of mass as:

\[E(C_i) = \int_{x \in C_{i}}\Big( f(x) - \lambda_{min}(C_i)\Big)dx\]

This exhibits monotonic behavior along branches of the dendrogram, so it cannot be used to compare stabilities of nested clusters. Instead, we can define the relative excess of mass:

\[E_R(C_i) = \int_{x \in C_{i}}\Big( \lambda_{max}(x, C_i) - \lambda_{min}(C_i)\Big)dx\]

where \(\lambda_{max}(C_i)\) is the density level at which \(C_i\) is split or disappears, and \(\lambda_{max}(x, C_i) = \min\{f(x), \lambda_{max}(C_i)\}\)

For an HDBSCAN hierarchy, where we have finite data set \(X\), cluster labels, and desnity thresholds associated with each hierarchical level, we can adapt the previous expression to define the stability of a cluster \(C_i\) as:

\[ S(C_i) = \sum_{x_j \in C_i} \Big( \lambda_{max}(x, C_i) - \lambda_{min}(C_i)\Big) = \sum_{x_j \in C_i} \Big(\frac{1}{\epsilon_{min}(x_j, C_i)} - \frac{1}{\epsilon_{max}(C_i)}\Big)\]

where \(\lambda_{min}(C_i)\) is the minimum density level at which \(C_i\) exists, \(\lambda_{max}(x_j, C_i)\) is the density level beyond which object \(x_j\) no longer belongs to cluster \(C_i\) and \(\epsilon_{max}(C_i)\) and \(\epsilon_{min}(x_j, C_i)\) are the corresponding values for the threshold \(\epsilon\).

Cluster Selection Algorithm

With the notion of cluster stability developed, we can have extraction of the most prominent clusters formulated as an optimization problem of maximizing cluster stabilities subject to constraints:

\[\underset{\delta_2, ..., \delta_{\kappa}}\max J = \sum^\kappa_{i=2}\delta_iS(C_i)\]

\[\text{subject to} \begin{cases} \delta_i \in \{0,1\}, i=2,...,\kappa \\ \sum_{j \in I_h} \delta_j = 1, \forall h \in L\\ \end{cases}\]

where \(\delta_i\) indicates if cluster \(i\) is in the flat solution, \(L = \{h | C_h \text{is a leaf cluster}\}\) is the set of indices of leaf clusters, and \(I_h = \{j | j \ne 1 \text{ and }C_j \text{ is an ascendant of } C_h\}\), the set of indices of all clusters on the path from \(C_h\) (included) to the root (excluded). The constraints prevent nested clusters on the same path from being selected at the same time.

Selecting the clusters requires bottom up processing of nodes excluding the root, starting with the leaves, deciding whether to keep node \(C_i\) or a best so far selection of clusters in its subtrees. In the process, the total stability, \(\hat{S}(C_i)\) is updated:

\[\hat{S}(C_i) = \begin{cases} S(C_i), &\text{if } C_i \text{ is a leaf node}\\ \max\{S(C_i), \hat{S}(C_{i_l}) + \hat{S}(C_{i_r}) &\text{if } C_i \text{ is an internal node}\\ \end{cases} \]

Steps

1. Initialize \(\delta_2 = ... \delta_\kappa = 1\) and for all leaf nodes, set \(\hat{S}(C_h) = S(C_h)\)
2. Starting from the deepest levels, do bottom-up (excluding root):
   • 2.1 If \(S(C_i) < \hat{S}(C_{i_r}) + \hat{S}(C_{i_l})\), set \(\hat{S}(C_i) = \hat{S}(C_{i_l}) + \hat{S}(C_{i_r})\) and set \(\delta_i=0\)
   • 2.2 Else, set \(\hat{S}(C_i) = S(C_i)\) and set \(\delta_{(.)} = 0\) for all clusters in \(C_i\)’s subtrees.

Implementation

Hierarchical clustering can typically be slow but thanks to parallelization, can be sped up via GPUs. Nvidia’s blog entry, Faster HDBSCAN Soft Clustering with RAPIDS cuML, shows comparisons with multiple methods on datasets of varying sizes.

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn.datasets as data
%matplotlib inline
sns.set_context('poster')
sns.set_style('white')
sns.set_color_codes()
plot_kwds = {'alpha' : 0.5, 's' : 80, 'linewidths':0}

moons, _ = data.make_moons(n_samples=50, noise=0.05)
blobs, _ = data.make_blobs(n_samples=50, centers=[(-0.75,2.25), (1.0, 2.0)], cluster_std=0.25)
test_data = np.vstack([moons, blobs])
plt.scatter(test_data.T[0], test_data.T[1], color='b', **plot_kwds)

import hdbscan

clusterer = hdbscan.HDBSCAN(min_cluster_size=5, gen_min_span_tree=True)
clusterer.fit(test_data)

HDBSCAN(gen_min_span_tree=True)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

clusterer.minimum_spanning_tree_.plot(edge_cmap='viridis',
                                      edge_alpha=0.6,
                                      node_size=80,
                                      edge_linewidth=2)

clusterer.single_linkage_tree_.plot(cmap='viridis', colorbar=True)

clusterer.condensed_tree_.plot()

clusterer.condensed_tree_.plot(select_clusters=True, selection_palette=sns.color_palette())

palette = sns.color_palette()
cluster_colors = [sns.desaturate(palette[col], sat)
                  if col >= 0 else (0.5, 0.5, 0.5) for col, sat in
                  zip(clusterer.labels_, clusterer.probabilities_)]
plt.scatter(test_data.T[0], test_data.T[1], c=cluster_colors, **plot_kwds)

Soft Clustering for HDBSCAN

Concept

• Hard clustering assigns a single cluster label or noise
• Soft clustering assigns a vector of probabilities

Instructional Implementation

import hdbscan
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import matplotlib as mpl

from scipy.spatial.distance import cdist

%matplotlib inline
sns.set_context('poster')
sns.set_style('white')
sns.set_color_codes()

plot_kwds={'alpha':0.25, 's':60, 'linewidths':0}
palette = sns.color_palette('deep', 12)

data = np.load('clusterable_data.npy')
fig = plt.figure()
ax = fig.add_subplot(111)
plt.scatter(data.T[0], data.T[1], **plot_kwds)
ax.set_xticks([])
ax.set_yticks([]);

len(data)

2309

# build a clustering object to fit the data
clusterer = hdbscan.HDBSCAN(min_cluster_size=15).fit(data)
clusterer.labels_

array([ 5,  5,  5, ..., -1, -1,  5], dtype=int64)

len(clusterer.labels_)

2309

clusterer.probabilities_

array([1.        , 0.85883269, 0.90828071, ..., 0.        , 0.        ,
       1.        ])

clusterer.probabilities_

array([1.        , 0.85883269, 0.90828071, ..., 0.        , 0.        ,
       1.        ])

# visualize clustering using scores as saturation
pal = sns.color_palette('deep', 8)
colors = [sns.desaturate(pal[col], sat) for col, sat in zip(clusterer.labels_,
                                                            clusterer.probabilities_)]
plt.scatter(data.T[0], data.T[1], c=colors, **plot_kwds)

Distance Based Membership

This takes advantage of features to examine clustering further, as pointed out in these docs. In particular, we use the CondensedTree object, obtained via `clusterer. which is a dendrogram, or family tree of various clusters.

• _raw_tree is the tree in the form of a numpy recarray allowing lookup based on field names where each row represents an edge between a parent and child clusters:
  • parent: id of parent cluster
  • child: id of child cluster
  • lambda_val: inverse distance, aka density at which edge forms
  • child_size: number of points in the child cluster

Outlier Based Membership

• these are density-based memberships using cluster persistence as a baseline of comparison
• GLOSH algorithm compares cluster persistence with how long a point stayed in the cluster
• we can compare cluster persistence with the merge height of the point with a fixed cluster in the dendrogram
  • merge height represents the disimilarity between two clusters
  • we can do that with multiple clusters to form a vector which can be normalized into memership scores

def max_lambda_val(cluster, tree):
    cluster_tree = tree[tree['child_size'] > 1]
    leaves = hdbscan.plots._recurse_leaf_dfs(cluster_tree, cluster)
    max_lambda = 0.0
    for leaf in leaves:
        max_lambda = max(max_lambda,
                         tree['lambda_val'][tree['parent'] == leaf].max())
    return max_lambda

def points_in_cluster(cluster, tree):
    leaves = hdbscan.plots._recurse_leaf_dfs(tree, cluster)
    return leaves

def merge_height(point, cluster, tree, point_dict):
    cluster_row = tree[tree['child'] == cluster]
    cluster_height = cluster_row['lambda_val'][0]
    if point in point_dict[cluster]:
        merge_row = tree[tree['child'] == float(point)][0]
        return merge_row['lambda_val']
        
    else:
        while point not in point_dict[cluster]:
            parent_row = tree[tree['child'] == cluster]
            cluster = parent_row['parent'].astype(np.float64)[0]
        for row in tree[tree['parent'] == cluster]:
            child_cluster = float(row['child'])
            if child_cluster == point:
                return row['lambda_val']
            if child_cluster in point_dict and point in point_dict[child_cluster]:
                return row['lambda_val']
            
def per_cluster_scores(point, cluster_ids, tree, max_lambda_dict, point_dict):
    result = {}
    point_row = tree[tree['child'] == point]
    point_cluster = float(point_row[0]['parent'])
    max_lambda = max_lambda_dict[point_cluster] + 1e-8 # avoid zero lambda vals in odd cases
    
    for c in cluster_ids:
        height = merge_height(point, c, tree, point_dict)
        result[c] = (max_lambda / (max_lambda - height))
    return result

def outlier_membership_vector(point, cluster_ids, tree,
                              max_lambda_dict, point_dict, softmax=True):
    if softmax:
        result = np.exp(np.array(list(per_cluster_scores(point,
                                                         cluster_ids,
                                                         tree,
                                                         max_lambda_dict,
                                                         point_dict
                                                         ).values())))
        result[~np.isfinite(result)] = np.finfo(np.double).max
    else:
        result = np.array(list(per_cluster_scores(point,
                                                  cluster_ids,
                                                  tree,
                                                  max_lambda_dict,
                                                  point_dict
                                                 ).values()))
        result /= result.sum()
    return result

cluster_ids = tree._select_clusters()
raw_tree = tree._raw_tree
all_possible_clusters = np.arange(data.shape[0], raw_tree['parent'].max() + 1).astype(np.float64)
max_lambda_dict = {c:max_lambda_val(c, raw_tree) for c in all_possible_clusters}
point_dict = {c:set(points_in_cluster(c, raw_tree)) for c in all_possible_clusters}
colors = np.empty((data.shape[0], 3))
for x in range(data.shape[0]):
    membership_vector = outlier_membership_vector(x, cluster_ids, raw_tree,
                                                  max_lambda_dict, point_dict, False)
    color = np.argmax(membership_vector)
    saturation = membership_vector[color]
    colors[x] = sns.desaturate(pal[color], saturation)
plt.scatter(data.T[0], data.T[1], c=colors, **plot_kwds);

Middle Way

To combine distance and membership approaches we can: - view membership vectors as probability distributions - combining them can be achieved via Bayes’ rule

Converting a Conditional Probability

The previous computation was the probability vector that a point is in each cluster, given that the point is in a cluster. To convert the conditional probability to a joint one, we need to: - multiply it by the probability that that there is a cluster to which the point belongs - this can be estimated from the merge height and comparing it with the max density - this should result in a number between 0 and 1

• plotting can show the $ $ cluster and the corresponding color

def prob_in_some_cluster(point, tree, cluster_ids, point_dict, max_lambda_dict):
    heights = []
    for cluster in cluster_ids:
        heights.append(merge_height(point, cluster, tree._raw_tree, point_dict))
    height = max(heights)
    nearest_cluster = cluster_ids[np.argmax(heights)]
    max_lambda = max_lambda_dict[nearest_cluster]
    return height / max_lambda

colors = np.empty((data.shape[0], 3))
for x in range(data.shape[0]):
    membership_vector = combined_membership_vector(x, data, tree, exemplar_dict, cluster_ids,
                                                   max_lambda_dict, point_dict, False)
    membership_vector *= prob_in_some_cluster(x, tree, cluster_ids, point_dict, max_lambda_dict)
    color = np.argmax(membership_vector)
    saturation = membership_vector[color]
    colors[x] = sns.desaturate(pal[color], saturation)
plt.scatter(data.T[0], data.T[1], c=colors, **plot_kwds);

C:\Users\Jonathan\AppData\Local\Temp\ipykernel_16588\1075745609.py:16: RuntimeWarning: divide by zero encountered in divide
  result = 1./dist_vector(point, exemplar_dict, data)

Soft Clustering in Practice

What are the general steps involved in performing a soft clustering analysis?

1. dimensionality reduction as a prerequisite to 2D visualization:
   • PCA
   • TSNE
   • UMAP
     
2. visualization of raw data
3. soft clustering via HDBSCAN
4. visualization of cluster membership
   • hard clustering color coding
   • soft clustering color saturation
     
5. quantitative analysis of cluster membership probabilities

from sklearn import datasets
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import hdbscan

# loading raw data, 8x8 gray scale of handwritten digits, plotting low dimensional projection
digits = datasets.load_digits()
data = digits.data
projection = TSNE().fit_transform(data)
plt.scatter(*projection.T, **plot_kwds)


# setting up a clustering object, fitting and clustering the data 
clusterer = hdbscan.HDBSCAN(min_cluster_size=10, prediction_data=True).fit(data)
color_palette = sns.color_palette('Paired', 12)
cluster_colors = [color_palette[x] if x >= 0
                  else (0.5, 0.5, 0.5)
                  for x in clusterer.labels_]
cluster_member_colors = [sns.desaturate(x, p) for x, p in
                         zip(cluster_colors, clusterer.probabilities_)]
plt.scatter(*projection.T, s=50, linewidth=0, c=cluster_member_colors, alpha=0.25)

Some of the data is noisy, so we can examine the noisy ones more through soft clustering.

soft_clusters = hdbscan.all_points_membership_vectors(clusterer)
color_palette = sns.color_palette('Paired', 12)
cluster_colors = [color_palette[np.argmax(x)]
                  for x in soft_clusters]
plt.scatter(*projection.T, s=50, linewidth=0, c=cluster_colors, alpha=0.25)

We can show uncertainty by coupling probability with desaturation of color. The higher the probability, the more saturated or pure the color. The lower the probability, the more gray it will be. Here we can see that desaturation is a harsh treatment with lots of gray, and that a lower limit may be more visually meaningful.

color_palette = sns.color_palette('Paired', 12)
cluster_colors = [sns.desaturate(color_palette[np.argmax(x)], np.max(x))
                  for x in soft_clusters]
plt.scatter(*projection.T, s=50, linewidth=0, c=cluster_colors, alpha=0.25)

One question to investigate is what points have high likelihoods for two clusters and low for the others. It’s worthwhile to note that the probabilities are joint ones, and that points have a probability of not being in a cluster.

def top_two_probs_diff(probs):
    sorted_probs = np.sort(probs)
    return sorted_probs[-1] - sorted_probs[-2]

# Compute the differences between the top two probabilities
diffs = np.array([top_two_probs_diff(x) for x in soft_clusters])
# Select out the indices that have a small difference, and a larger total probability, extract from tuple form
mixed_points = np.where((diffs < 0.001)  & (np.sum(soft_clusters, axis=1) > 0.5))[0]

colors = [(0.75, 0.1, 0.1) if x in mixed_points
          else (0.5, 0.5, 0.5) for x in range(data.shape[0])]
plt.scatter(*projection.T, s=50, linewidth=0, c=colors, alpha=0.5)

fig = plt.figure()
for i, image in enumerate(digits.images[mixed_points][:16]):
    ax = fig.add_subplot(4,4,i+1)
    ax.imshow(image)
plt.tight_layout()

plt.hist(np.sum(soft_clusters, axis=1))

(array([ 36.,  47.,  92., 171., 332., 229., 181., 297., 153., 259.]),
 array([0.58730357, 0.62857321, 0.66984285, 0.7111125 , 0.75238214,
        0.79365178, 0.83492143, 0.87619107, 0.91746071, 0.95873036,
        1.        ]),
 <BarContainer object of 10 artists>)

K-means Clustering

Description

Given a set of observations \((x_1, x_2, ..., x_n)\) in a \(d\)-dimensional vector space, k-means clustering partitions the \(n\) observations into \(k\) sets, \(\textbf{S} = \{S_1, S_2, ... S_k\}\). The goal is to minimize the variance within each cluster:

\[\underset{\textbf{S}}{\arg\min} \sum^k_{i=1} \sum_{x\in S_i} ||x - \mu_i||^2 = \underset{\textbf{S}}{\arg\min}\sum^k_{i=1}|S_i|\text{Var}S_i\]

where \(\mu_i\) is the mean or centroid of the points in \(S_i\).

Steps

Initialization

First choose \(k\) and initialize randomly chosen \(k\) centroids at \(t = 0\), the initial step.

Assignment

Assign each observation to the cluster with the nearest mean, in terms of Euclidean distance. Each set is thus:

\[S^{(t)}_i = \big\{x_p : ||x_p - m_i^{(t)}||^2 \le ||x_p - m_j^{(t)}||^2 \forall j, 1 \le j \le k \big\}\]

Update

After assignment, update to find the new means for the next iteration, \(t + 1\):

\[ m_i^{t + 1} = \frac{1}{|S_i^{(t)}|}\sum_{x_j \in S_i^{(t)}}x_j\]

The number of computations scales with the number of samples \(n\), the number of centers \(k\), the number of dimensions \(d\) and the number of iterations \(i\), for \(O(nkdi)\)

Creating the Dataset

# import statements
from sklearn.datasets import make_blobs
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set('notebook')
# Generate and plot data with predetermined model and classes
X, y = make_blobs(n_samples=200, n_features=2, centers=4,
                  cluster_std=1.6, random_state=50)

plt.scatter(X[:,0], X[:,1], c=y, cmap='viridis')
plt.xlim(-15,15)
plt.ylim(-15,15)
plt.show()

Y = np.zeros((200, 2))
Y[:,0] = y
Y

array([[2., 0.],
       [0., 0.],
       [3., 0.],
       [3., 0.],
       [1., 0.],
       [0., 0.],
       [2., 0.],
       [0., 0.],
       [0., 0.],
       [1., 0.],
       [1., 0.],
       [2., 0.],
       [1., 0.],
       [1., 0.],
       [1., 0.],
       [3., 0.],
       [2., 0.],
       [2., 0.],
       [1., 0.],
       [1., 0.],
       [1., 0.],
       [3., 0.],
       [2., 0.],
       [2., 0.],
       [2., 0.],
       [3., 0.],
       [2., 0.],
       [1., 0.],
       [3., 0.],
       [2., 0.],
       [2., 0.],
       [1., 0.],
       [1., 0.],
       [0., 0.],
       [0., 0.],
       [1., 0.],
       [1., 0.],
       [1., 0.],
       [3., 0.],
       [3., 0.],
       [2., 0.],
       [1., 0.],
       [3., 0.],
       [3., 0.],
       [0., 0.],
       [0., 0.],
       [2., 0.],
       [3., 0.],
       [3., 0.],
       [2., 0.],
       [0., 0.],
       [1., 0.],
       [3., 0.],
       [2., 0.],
       [3., 0.],
       [3., 0.],
       [0., 0.],
       [1., 0.],
       [1., 0.],
       [3., 0.],
       [0., 0.],
       [0., 0.],
       [2., 0.],
       [0., 0.],
       [0., 0.],
       [2., 0.],
       [0., 0.],
       [3., 0.],
       [2., 0.],
       [2., 0.],
       [0., 0.],
       [1., 0.],
       [2., 0.],
       [0., 0.],
       [3., 0.],
       [3., 0.],
       [1., 0.],
       [3., 0.],
       [3., 0.],
       [0., 0.],
       [1., 0.],
       [3., 0.],
       [1., 0.],
       [1., 0.],
       [1., 0.],
       [1., 0.],
       [2., 0.],
       [0., 0.],
       [2., 0.],
       [1., 0.],
       [0., 0.],
       [2., 0.],
       [0., 0.],
       [3., 0.],
       [3., 0.],
       [1., 0.],
       [0., 0.],
       [3., 0.],
       [0., 0.],
       [3., 0.],
       [1., 0.],
       [0., 0.],
       [3., 0.],
       [2., 0.],
       [3., 0.],
       [0., 0.],
       [3., 0.],
       [1., 0.],
       [1., 0.],
       [2., 0.],
       [0., 0.],
       [1., 0.],
       [2., 0.],
       [2., 0.],
       [1., 0.],
       [3., 0.],
       [2., 0.],
       [2., 0.],
       [3., 0.],
       [0., 0.],
       [2., 0.],
       [2., 0.],
       [2., 0.],
       [2., 0.],
       [0., 0.],
       [2., 0.],
       [2., 0.],
       [3., 0.],
       [2., 0.],
       [0., 0.],
       [1., 0.],
       [0., 0.],
       [0., 0.],
       [3., 0.],
       [2., 0.],
       [3., 0.],
       [0., 0.],
       [2., 0.],
       [2., 0.],
       [1., 0.],
       [2., 0.],
       [3., 0.],
       [2., 0.],
       [1., 0.],
       [1., 0.],
       [2., 0.],
       [0., 0.],
       [3., 0.],
       [0., 0.],
       [1., 0.],
       [2., 0.],
       [1., 0.],
       [3., 0.],
       [0., 0.],
       [0., 0.],
       [3., 0.],
       [2., 0.],
       [0., 0.],
       [0., 0.],
       [0., 0.],
       [3., 0.],
       [3., 0.],
       [1., 0.],
       [0., 0.],
       [1., 0.],
       [1., 0.],
       [1., 0.],
       [3., 0.],
       [0., 0.],
       [0., 0.],
       [0., 0.],
       [0., 0.],
       [1., 0.],
       [3., 0.],
       [0., 0.],
       [2., 0.],
       [1., 0.],
       [0., 0.],
       [2., 0.],
       [3., 0.],
       [0., 0.],
       [3., 0.],
       [3., 0.],
       [2., 0.],
       [1., 0.],
       [0., 0.],
       [1., 0.],
       [2., 0.],
       [3., 0.],
       [0., 0.],
       [3., 0.],
       [3., 0.],
       [1., 0.],
       [3., 0.],
       [0., 0.],
       [2., 0.],
       [1., 0.],
       [3., 0.],
       [1., 0.],
       [2., 0.]])

from sklearn.cluster import KMeans

n_clusters=4

# Fit and predict clusters
kmeans = KMeans(n_clusters=n_clusters)
kmeans.fit(X)
print(kmeans.cluster_centers_)
y_km = kmeans.fit_predict(X)

# Plot predictions, one class at a time
for i, c in enumerate(['red', 'black', 'blue', 'cyan']):
    plt.scatter(X[y_km ==i,0], X[y_km == i,1], s=100, c=c)
plt.show()

[[-5.56465793 -2.34988939]
 [-1.92101646  5.21673484]
 [ 0.05161133 -5.35489826]
 [-2.40167949 10.17352695]]

Hierarchical Clustering

Hierarchical clustering involves starting with treating each observation as a set, and then at the following step, creating a new cluster from the two “nearest” clusters, according to the defined distance metric. There are multiple ways the idea of a distance metric can extend to clusters. Some common choices include the average, closest members, or farthest nembers.

# import hierarchical clustering libraries
from sklearn.cluster import AgglomerativeClustering
from matplotlib import pyplot as plt
import seaborn as sns
from scipy.cluster.hierarchy import dendrogram, linkage, to_tree
%load_ext autoreload
%autoreload 2
sns.set()
data = [[i] for i in [-30, 4, 1, 2, 5, 6, 10, 50, 75, 100]]

def llf(id):
    return str(id)

Z = linkage(data, 'ward')
fig = plt.figure(figsize=(20, 10))
dn_truncated = dendrogram(Z, orientation='right', p=4, truncate_mode='lastp', leaf_label_func=llf)

dn = dendrogram(Z, orientation='right', p=4, truncate_mode=None, leaf_label_func=llf)

dn_truncated

{'icoord': [[5.0, 5.0, 15.0, 15.0],
  [25.0, 25.0, 35.0, 35.0],
  [10.0, 10.0, 30.0, 30.0]],
 'dcoord': [[0.0, 43.30127018922193, 43.30127018922193, 0.0],
  [0.0, 45.389321169086415, 45.389321169086415, 0.0],
  [43.30127018922193,
   154.2898015331631,
   154.2898015331631,
   45.389321169086415]],
 'ivl': ['9', '15', '0', '14'],
 'leaves': [9, 15, 0, 14],
 'color_list': ['C1', 'C2', 'C0'],
 'leaves_color_list': ['C1', 'C1', 'C2', 'C2'],
 'traversal': [18, 16, 9, 15, 17, 0, 14]}

dn['icoord']

[[15.0, 15.0, 25.0, 25.0],
 [5.0, 5.0, 20.0, 20.0],
 [55.0, 55.0, 65.0, 65.0],
 [85.0, 85.0, 95.0, 95.0],
 [75.0, 75.0, 90.0, 90.0],
 [60.0, 60.0, 82.5, 82.5],
 [45.0, 45.0, 71.25, 71.25],
 [35.0, 35.0, 58.125, 58.125],
 [12.5, 12.5, 46.5625, 46.5625]]

dn['dcoord']

[[0.0, 25.0, 25.0, 0.0],
 [0.0, 43.30127018922193, 43.30127018922193, 25.0],
 [0.0, 1.0, 1.0, 0.0],
 [0.0, 1.0, 1.0, 0.0],
 [0.0, 1.7320508075688772, 1.7320508075688772, 1.0],
 [1.0, 5.422176684690383, 5.422176684690383, 1.7320508075688772],
 [0.0, 8.262364471909155, 8.262364471909155, 5.422176684690383],
 [0.0, 45.389321169086415, 45.389321169086415, 8.262364471909155],
 [43.30127018922193, 154.2898015331631, 154.2898015331631, 45.389321169086415]]

dn_truncated['icoord']

[[5.0, 5.0, 15.0, 15.0], [25.0, 25.0, 35.0, 35.0], [10.0, 10.0, 30.0, 30.0]]

dn_truncated['dcoord']

[[0.0, 43.30127018922193, 43.30127018922193, 0.0],
 [0.0, 45.389321169086415, 45.389321169086415, 0.0],
 [43.30127018922193, 154.2898015331631, 154.2898015331631, 45.389321169086415]]

k = 4
for i in zip(dn['icoord'][k], dn['dcoord'][k]):
    print(i)

(75.0, 0.0)
(75.0, 1.7320508075688772)
(90.0, 1.7320508075688772)
(90.0, 1.0)

Linkage Matrix

Full documentation for Scipy’s implementation is found here. The original observations are treated as \(n\) clusters. There are \(n-1\) new clusters created beyond these total, with indices ranging from \(n+1\) to \(2n-1\), ordered by the inter-cluster distance. Each of thse new clusters is made from two other clusters. The \(n-1 \times 4\) linkage matrix \(Z\) encodes information about the new clusters in the following manner: Each row \(i\) (keeping in mind zero indexing) describes the \((i+1)\)th cluster by listing the indices of the two source clusters, the inter-cluster distance between these two clusters, and finally, the number of members of the new cluster.

Z

array([[  2.        ,   3.        ,   1.        ,   2.        ],
       [  1.        ,   4.        ,   1.        ,   2.        ],
       [  5.        ,  11.        ,   1.73205081,   3.        ],
       [ 10.        ,  12.        ,   5.42217668,   5.        ],
       [  6.        ,  13.        ,   8.26236447,   6.        ],
       [  7.        ,   8.        ,  25.        ,   2.        ],
       [  9.        ,  15.        ,  43.30127019,   3.        ],
       [  0.        ,  14.        ,  45.38932117,   7.        ],
       [ 16.        ,  17.        , 154.28980153,  10.        ]])

plt.show()

Converting the Linkage Matrix to Tree

to_tree(Z)
rootnode, nodelist = to_tree(Z, rd=True)
rootnode

<__main__.ClusterNode at 0x29427b86370>

rootnode.get_right().pre_order()

[0, 6, 2, 3, 5, 1, 4]

len(nodelist)

19

nodelist[4].pre_order()

[4]

Visualization with Plotly

p =6

np.array(data)

array([[-30],
       [  4],
       [  1],
       [  2],
       [  5],
       [  6],
       [ 10],
       [ 50],
       [ 75],
       [100]])

n = len(data)

len(range(2*n-p + 1, 2*n))

5

traversal = dn_truncated['traversal']

traversal.reverse()

len(traversal)

7

traversal[p:]

[14]

# import plotly.figure_factory as ff
import numpy as np
np.random.seed(1)
n = len(data)
p = 4
fig = create_dendrogram(np.array(data), hovertext=traversal[p:],
                        labels=data, orientation='left', p=p, truncate_mode='lastp', leaf_label_func=llf)
fig.update_layout(width=800, height=500)
fig.show()

nodelist[13].pre_order()

[2, 3, 5, 1, 4]

dn_truncated

{'icoord': [[15.0, 15.0, 25.0, 25.0],
  [5.0, 5.0, 20.0, 20.0],
  [45.0, 45.0, 55.0, 55.0],
  [35.0, 35.0, 50.0, 50.0],
  [12.5, 12.5, 42.5, 42.5]],
 'dcoord': [[0.0, 25.0, 25.0, 0.0],
  [0.0, 43.30127018922193, 43.30127018922193, 25.0],
  [0.0, 8.262364471909155, 8.262364471909155, 0.0],
  [0.0, 45.389321169086415, 45.389321169086415, 8.262364471909155],
  [43.30127018922193,
   154.2898015331631,
   154.2898015331631,
   45.389321169086415]],
 'ivl': ['9', '7', '8', '0', '6', '13'],
 'leaves': [9, 7, 8, 0, 6, 13],
 'color_list': ['C1', 'C1', 'C2', 'C2', 'C0'],
 'leaves_color_list': ['C1', 'C1', 'C1', 'C2', 'C2', 'C2']}

Z

array([[  2.        ,   3.        ,   1.        ,   2.        ],
       [  1.        ,   4.        ,   1.        ,   2.        ],
       [  5.        ,  11.        ,   1.73205081,   3.        ],
       [ 10.        ,  12.        ,   5.42217668,   5.        ],
       [  6.        ,  13.        ,   8.26236447,   6.        ],
       [  7.        ,   8.        ,  25.        ,   2.        ],
       [  9.        ,  15.        ,  43.30127019,   3.        ],
       [  0.        ,  14.        ,  45.38932117,   7.        ],
       [ 16.        ,  17.        , 154.28980153,  10.        ]])

dn

{'icoord': [[15.0, 15.0, 25.0, 25.0],
  [5.0, 5.0, 20.0, 20.0],
  [55.0, 55.0, 65.0, 65.0],
  [85.0, 85.0, 95.0, 95.0],
  [75.0, 75.0, 90.0, 90.0],
  [60.0, 60.0, 82.5, 82.5],
  [45.0, 45.0, 71.25, 71.25],
  [35.0, 35.0, 58.125, 58.125],
  [12.5, 12.5, 46.5625, 46.5625]],
 'dcoord': [[0.0, 25.0, 25.0, 0.0],
  [0.0, 43.30127018922193, 43.30127018922193, 25.0],
  [0.0, 1.0, 1.0, 0.0],
  [0.0, 1.0, 1.0, 0.0],
  [0.0, 1.7320508075688772, 1.7320508075688772, 1.0],
  [1.0, 5.422176684690383, 5.422176684690383, 1.7320508075688772],
  [0.0, 8.262364471909155, 8.262364471909155, 5.422176684690383],
  [0.0, 45.389321169086415, 45.389321169086415, 8.262364471909155],
  [43.30127018922193,
   154.2898015331631,
   154.2898015331631,
   45.389321169086415]],
 'ivl': ['9', '7', '8', '0', '6', '2', '3', '5', '1', '4'],
 'leaves': [9, 7, 8, 0, 6, 2, 3, 5, 1, 4],
 'color_list': ['C1', 'C1', 'C2', 'C2', 'C2', 'C2', 'C2', 'C2', 'C0'],
 'leaves_color_list': ['C1',
  'C1',
  'C1',
  'C2',
  'C2',
  'C2',
  'C2',
  'C2',
  'C2',
  'C2'],
 'traversal': [18,
  16,
  9,
  15,
  7,
  8,
  17,
  0,
  14,
  6,
  13,
  10,
  2,
  3,
  12,
  5,
  11,
  1,
  4]}

dn_truncated

{'icoord': [[15.0, 15.0, 25.0, 25.0],
  [5.0, 5.0, 20.0, 20.0],
  [45.0, 45.0, 55.0, 55.0],
  [35.0, 35.0, 50.0, 50.0],
  [12.5, 12.5, 42.5, 42.5]],
 'dcoord': [[0.0, 25.0, 25.0, 0.0],
  [0.0, 43.30127018922193, 43.30127018922193, 25.0],
  [0.0, 8.262364471909155, 8.262364471909155, 0.0],
  [0.0, 45.389321169086415, 45.389321169086415, 8.262364471909155],
  [43.30127018922193,
   154.2898015331631,
   154.2898015331631,
   45.389321169086415]],
 'ivl': ['9', '7', '8', '0', '6', '(5)'],
 'leaves': [9, 7, 8, 0, 6, 13],
 'color_list': ['C1', 'C1', 'C2', 'C2', 'C0'],
 'leaves_color_list': ['C1', 'C1', 'C1', 'C2', 'C2', 'C2']}

i = 0
print(dn_truncated['dcoord'][i])
print(dn_truncated['icoord'][i])

[0.0, 25.0, 25.0, 0.0]
[15.0, 15.0, 25.0, 25.0]

dn_truncated['dcoord']

[[0.0, 25.0, 25.0, 0.0],
 [0.0, 43.30127018922193, 43.30127018922193, 25.0],
 [0.0, 8.262364471909155, 8.262364471909155, 0.0],
 [0.0, 45.389321169086415, 45.389321169086415, 8.262364471909155],
 [43.30127018922193, 154.2898015331631, 154.2898015331631, 45.389321169086415]]

dn_truncated['icoord']

[[15.0, 15.0, 25.0, 25.0],
 [5.0, 5.0, 20.0, 20.0],
 [45.0, 45.0, 55.0, 55.0],
 [35.0, 35.0, 50.0, 50.0],
 [12.5, 12.5, 42.5, 42.5]]

nodelist[].pre_order()

[4]

nodelist[9].id

data[2]

[1]

#from plotly.figure_factory import create_dendrogram
import numpy as np
import pandas as pd
Index= ['A','B','C','D','E','F','G','H','I','J']
df = pd.DataFrame(abs(np.random.randn(10, 10)), index=Index)
fig = create_dendrogram(df, labels=Index, orientation='left', hovertext=list(range(19)))
fig.show()

Implementing Cluster Order in Scipy

Modification of _dendrogram_calculate_info function

• _dendrogram_calculate_info now has a list to track the cluster indices during construction of the dendrogram

Implementing Dendrogram Cuts in Plotly

Modification of create_dendrogram function and _Dendrogram class

To implement cuts, the basic modification to be made is including truncation parameters such as truncate_mode and p to all wrappers around the scipy’s dendrogram function. The following changes are made:

• create_dendrogram now includes p and truncate_mode as parameters

• the plotly class _Dendrogram now has both p and truncate_mode as object parameters

• the get_dendrogram_traces method now passes on these parameters to its call to scipy’s dendogram function

Creating Hovertext Robust to Truncation

• In Dendrograms with hovertext ranging a few hundred nodes, truncation is necessary

• The Z matrix contains contains \(n-1\) rows each representing new clusters, where \(n\) is then number of data points.

  • higher rows signify larger clusters, ordered by the distance metric of their child clusters.
  • the row indices also correspond to the ClusterNode instances in the node list returned by scipy.cluster.hierarchy.to_tree
  • default indexing will have the root cluster as \(n-2\)

• Truncation when truncate_mode is set to ‘lastp’ is governed by parameter p

  • p becomes the number of leaves to be displayed
  • there will be a total of \(p - 1\) new clusters displayed
  • the node indices will then be in range [2n-2-p, 2n-2], or range(2*n - 2 - p, 2*n - 1)

• If indices are shown

• For UI considerations, a shift in p Simpler indices, \(j\) can be transformations of default indices \(i\) for UI considerations can be given by: \[j = 2n-2 - i\]

Visualization Techniques

t-Stochastic Neighbor Embedding (TSNE)

This is a technique that maps a higher dimensional space to a lower dimensional space while preserving probabilities of choosing analogous neighboring points between the two difference spaces.

The original paper can be read here as well as a post about effective usage.

Uniform Manifold Approximation (UMAP)

This is another method of dimensional reduction but works well for non-linear mappings.

Implementation with Plotly

Here is an example of using Plotly to visualize data sets using both techniques, straight from their documentation

import plotly.express as px

df = px.data.iris()
features = ["sepal_width", "sepal_length", "petal_width", "petal_length"]
fig = px.scatter_matrix(df, dimensions=features, color="species")
fig.show()

from sklearn.manifold import TSNE, trustworthiness
import plotly.express as px

df = px.data.iris()

features = df.loc[:, :'petal_width']

tsne = TSNE(n_components=2, random_state=0)
projections = tsne.fit_transform(features)

fig = px.scatter(
    projections, x=0, y=1,
    color=df.species, labels={'color': 'species'}
)
fig.show()

score = trustworthiness(features.to_numpy(), projections, n_neighbors=5)
print(f'Trustworthiness score is: {score}')

Trustworthiness score is: 0.9880657276995305

from umap import UMAP
import plotly.express as px

df = px.data.iris()

features = df.loc[:, :'petal_width']
umap_2d = UMAP(random_state=0, n_components=2)
umap_3d = UMAP(random_state=0, n_components=3, n_neighbors=5)
proj_2d = umap_2d.fit_transform(features)
fig_2d = px.scatter(proj_2d, x=0, y=1,
                 color=df.species, labels={'color': 'species'}
                )
fig_2d.show()

C:\Users\Jonathan\anaconda3\envs\py313\Lib\site-packages\umap\umap_.py:1952: UserWarning:

n_jobs value 1 overridden to 1 by setting random_state. Use no seed for parallelism.

score = trustworthiness(features.to_numpy(), projections, n_neighbors=5)
print(f'Trustworthiness score is: {score}')

Trustworthiness score is: 0.9880657276995305